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Preface 


The  work  reported  herein  was  conducted  as  part  of  the  Water  Quality  Research  Program 
(WQRP),  Work  Unit  Number  33251.  The  WQRP  is  sponsored  by  Headquarters,  U.S.  Army  Corps 
of  Engineers  (HQUSACE),  and  is  assigned  to  the  U.S.  Army  Engineer  Research  and  Development 
Center  (ERDC)  under  the  purview  of  the  Environmental  Laboratory  (EL),  Vicksburg,  MS. 
Dr.  John  W.  Barko  was  Program  Manager  for  the  WQRP,  and  Mr.  Robert  C.  Gunkel,  Jr.,  was 
Assistant  Manager.  Program  Monitor  during  this  study  was  Mr.  Frederick  B.  Juhle,  HQUSACE. 

This  report  presents  the  theoretical  basis  for  the  extension  of  CE-QUAL-W2  to  model  entire 
river  basins  including  steeply  sloping  rivers,  reservoirs,  and  estuaries. 

The  study  was  conducted  and  the  report  written  by  Dr.  Scott  A.  Wells  of  the  Department  of 
Civil  Engineering,  Portland  State  University,  Portland,  OR,  and  Mr.  Thomas  M.  Cole  of  the  Water 
Quality  and  Contaminant  Modeling  Branch  (WQCMB),  Environmental  Processes  and  Effects 
Division  (EPED),  EL.  Technical  review  by  Dr.  Mark  S.  Dortch,  Chief,  WQCMB,  and  Ms.  Lillian 
T.  Schneider,  WQCMB,  is  gratefully  acknowledged. 

This  investigation  was  conducted  under  the  general  supervision  of  Dr.  John  Keeley,  Acting 
Director,  EL,  and  Dr.  Richard  E.  Price,  Chief,  EPED,  and  under  the  direct  supervision  of 
Dr.  Dortch. 

At  the  time  of  publication  of  this  report.  Dr.  James  R.  Houston  was  Director  of  ERDC,  and 
COL  James  S.  Weller,  EN,  was  Commander. 

This  report  should  be  cited  as  follows: 

Wells,  S.  A.,  and  Cole,  T.  M.  (2000).  “Theoretical  basis  for  the  CE-QUAL-W2 
river  basin  model,”  ERDC/EL  TR-00-7,  U.S.  Army  Engineer  Research  and 
Development  Center,  Vicksburg,  MS. 
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Introduction 


CE-QUAL-W2  (W2)  is  a  two-dimensional  (2-D)  water  quality  and  hydrodynamic  code  supported 
by  the  USACE  Waterways  Experiments  Station  (Cole  and  Buchak,  1995).  The  model  has  been 
widely  applied  to  surface  water  systems  such  as  lakes,  reservoirs,  and  estuaries.  W2  predicts  water 
levels,  horizontal/vertical  velocities,  temperature,  and  21  other  water  quality  parameters.  A  typical 
grid  for  this  model  is  shown  in  Figure  1  where  the  vertical  axis  is  aligned  with  gravity. 


Two-dimensional  hydrodynamics 


Figure  1.  Typical  W2  grid. 


The  primary  objective  of  this  research  is  to  integrate  a  riverine  model  into  the  existing  W2  code  that 

would  provide  the  capability  for  modeling  entire  watersheds.  This  task  was  accomplished  by  the 

following  steps: 

•  Formal  derivation  of  governing  equations  and  solution  algorithm  with  general  channel 
slope 

•  Detailed  analysis  of  algorithm  for  linking  branches  and  smooth  implementation  of 
boundary  conditions  between  branches 

•  Algorithm  development  and  changes  to  basic  model  code  (including  branch  definitions 
with  slope,  slope  correction  to  solution  algorithm,  transfer  of  momentum  between 
internal  branches) 

These  tasks  were  performed  with  the  following  constraints  and  initiatives: 

•  Utilize  the  same  solution  algorithms  as  the  existing  code  for  hydrodynamics  and  water 
quality  for  the  riverine  system 
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•  Allow  momentum  transfer  between  adjacent  branches  for  internal  head  boundary 
conditions 

•  Refine  the  turbulence  closure  hypothesis  for  riverine  sections 


River  Basin  Modei  Development  Rationale 

W2  has  been  in  use  for  the  last  several  decades  as  a  tool  for  water  quality  managers  to  assess  the 
impacts  of  management  strategies  on  river,  reservoir,  lake,  and  estuarine  systems.  A  predominant 
feature  of  the  model  is  its  ability  to  compute  the  2-D  velocity  field  for  narrow  systems  that  stratify. 
In  contrast  to  many  reservoir  models  that  are  zero-dimensional  hydrodynamic  models,  an 
understanding  of  the  fluid  mechanical  transport  can  be  as  important  as  the  reaction  kinetics  in  the 
water  column. 

One  limitation  of  W2  is  its  inability  to  model  steeply  sloping  river  stretches  and  hence  an  entire 
water  basin.  Models,  such  as  WQRSS,  HEC-5Q,  and  HSPF  have  been  developed  for  water  basin 
modeling  but  have  serious  limitations.  The  most  serious  limitation  is  that  the  HEC-5Q  (similar  to 
WQRSS)  and  HSPF  models  incorporate  a  one-dimensional  (1-D)  longitudinal  river  model  with  a 
one-dimensional  vertical  reservoir  model  (1-D  in  water  quality  and  0-D  in  hydrodynamics).  The 
modeler  must  choose  the  location  of  the  transition  from  1-D  longitudinal  to  1-D  vertical.  Besides 
the  problem  of  not  solving  for  the  velocity  field  in  the  stratified,  reservoir  system,  any  point  source 
inputs  to  the  reservoir  section  are  spread  over  the  entire  longitudinal  distribution  of  the  reservoir 
cell.  This  has  created  problems  in  the  following  two  WQRSS  water  quality  modeling  studies. 

Wahiawa  Reservoir 

Wahiawa  Reservoir  (a  narrow,  5  mile  long  reservoir  with  100  ft  depth  at  the  dam).  The  HEC 
WQRSS  model  was  initially  applied  to  this  two-fork  reservoir  system.  The  system  is  shown  in 
Figure  2. 
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Figure  2.  Wahiawa  Reservoir,  Oahu,  Hawaii. 


The  WQRRS  model  schematization  is  contrasted  to  the  W2  schematization  for  Wahiawa  Reservoir 
in  Figure  3.  The  initial  reservoir  study  using  WQRRS  produced  poor  results  even  after  expending 
large  resources  to  “make”  the  model  work.  The  modeling  effort  did  not  provide  a  management  tool 
for  water  quality  managers  because  of  gross  errors  in  setting  up  the  model,  i.e.,  combining  the  two 
forks  and  spreading  the  discharge  from  the  wastewater  treatment  plant  throughout  the  length  of  the 
reservoir. 

Tualatin  River 

The  Tualatin  River,  located  in  Oregon,  is  a  32-mile  long,  narrow,  stratified  system,  with  pools  25- 
30  ft  deep.  The  WQRSS  model  was  applied  to  this  system  incorrectly  because  the  modelers 
decided  to  break  the  system  from  a  river  to  a  reservoir  at  the  location  of  a  wastewater  treatment 
plant  discharge.  Hence,  a  large  section  of  the  Tualatin  that  stratified  was  modeled  as  completely 
mixed  because  the  modelers  knew  it  would  be  inappropriate  to  spread  a  point  source  over  32  miles 
if  this  section  was  chosen  as  a  stratified  system.  A  later  application  of  W2  (Berger  and  Wells, 
1995)  correctly  represented  the  physics  of  the  system. 

In  these  two  cases,  the  application  of  WQRSS  had  serious  limitations  for  the  reservoir  section.  W2 
was  subsequently  applied  to  these  cases  and  was  able  to  be  used  effectively  because  of  its  2-D 
hydrodynamics  and  water  quality. 

Other  hydraulic  and  water  quality  models  in  common  use  for  unsteady  flow  include  the  1-D 
dynamic  EPA  model  DYNHYD  (Ambrose,  et  al.,  1 988),  used  together  with  the  multidimensional 
water  quality  model  WASP.  WASP  relies  on  DYNHYD  for  the  1-D  hydrodynamics.  If  WASP  is 
used  in  a  multi-dimensional  schematization,  the  modeler  must  supply  dispersion  coefficients  to 
allow  transport  in  the  vertical  or  lateral  directions.  In  addition,  the  Corps  model,  CE-QUAL-RIV1 
(Environmental  Laboratory,  1995),  is  a  one-dimensional  dynamic  flow  and  water  quality  model 
used  for  one-dimensional  river  or  stream  sections.  Each  of  these  models  does  not  have  the  ability  to 
characterize  adequately  the  hydraulics  or  water  quality  of  deeper  reservoir  systems  or  deep  river 
pools  that  stratify. 
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W2,  even  though  able  to  handle  narrow  systems  that  stratify,  is  not  well-suited  for  one-dimensional 
river  channels.  In  the  development  of  W2,  vertical  accelerations  were  considered  negligible 
compared  to  gravity  forces.  This  assumption  lead  to  the  hydrostatic  pressure  approximation  for  the 
z-momentum  equation.  In  sloping  channels,  this  assumption  is  not  always  valid  because  the  vertical 
accelerations  cannot  be  neglected  if  the  x  and  z  axes  are  aligned  with  an  elevation  datum  and 
gravity,  respectively.  In  addition,  the  current  W2  algorithm  does  not  allow  the  upstream  bed 
elevation  to  be  above  the  downstream  water  surface  elevation.  If  one  wanted  to  use  the  existing 
model  for  sloping  channels,  one  would  have  to  break  the  sloping  section  into  multiple  small 
branches.  Because  water  basin  modeling  is  becoming  more  and  more  essential  for  water  quality 
managers,  providing  the  capability  for  W2  to  be  used  as  a  complete  tool  for  water  basin  modeling  is 
an  essential  step  in  upgrading  the  state-of-the-art  in  modeling  river  basins. 
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WQRSS  Model  Schematic 


Inflows  from  wastewater  treatment 
plants,  North  and  South  Fork  of  K. 
streams 


CE-QUAL-W2  Schematic 


Figure  3.  Comparison  of  WQRSS  and  W2  schematization  for  Wahiawa  Reservoir. 


10 


Development  Approach 

There  are  many  approaches  that  could  be  implemented  within  W2  for  riverine  branches.  By 
choosing  a  theoretical  basis  for  the  riverine  branches  that  uses  the  existing  2-D  computational 
scheme  for  hydraulics  and  water  quality,  the  following  benefits  accrue: 

•  code  updates  in  the  computational  scheme  will  affect  the  entire  model  rather  than  just 
one  of  the  computational  schemes  for  either  the  riverine  or  the  reservoir  sections  leading 
to  easier  code  maintenance 

•  no  changes  would  be  made  to  the  temperature  or  water  quality  solution  algorithms 

•  by  using  the  two-dimensional  framework,  the  riverine  branches  would  also  have  the 
ability  to  predict  the  velocity  and  water  quality  field  in  two  dimensions.  This  has 
advantages  in  modeling  the  following  processes:  sediment  deposition  and  scour, 
particulate  (algae,  detritus,  suspended  solids)  sedimentation,  and  sediment  flux  processes. 

•  since  the  entire  watershed  model  has  the  same  theoretical  basis,  setting  up  branches  and 
interfacing  branches  involves  the  same  process  whether  for  reservoir  or  riverine  sections, 
thus  making  code  maintenance  and  model  set-up  easier. 

The  theoretical  approach  will  be  to  allow  each  branch  segment  to  have  a  channel  slope.  The 
governing  equations  will  then  be  re-derived  assuming  that  the  gravity  force  in  the  x  and  z- 
momentum  equations  is  adjusted  by  the  channel  slope.  This  is  shown  schematically  in  Figure  3. 


Figure  4.  Schematic  of  river-reservoir  linkage  where  a  is  the  slope  of  the  channel  bottom. 

The  turbulence  closure  hypothesis  used  in  the  reservoir  sections  may  not  be  completely  adequate  in 
the  riverine  branches.  Various  turbulence  closure  hypotheses  will  be  explored. 
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At  the  interface  between  the  two  branches  (a  riverine  and  a  reservoir  branch),  an  internal  head 
boundary  condition  would  be  specified.  In  the  present  model  configuration,  momentum  is  not 
transferred  between  branches  if  an  internal  head  boundary  condition  is  specified.  The  existing  code 
would  then  be  altered  to  allow  momentum  transfer  between  adjacent  branches  based  on  the  angle 
between  the  branches. 

As  a  test  case,  the  Tualatin  River  upper  reaches  will  be  used  with  the  new  algorithm.  Figure  5 
shows  the  river  bottom  elevations  as  a  function  of  river  mile.  Prior  modeling  of  the  pool  section 
was  reported  in  Berger  and  Wells  (1995). 


Diversion 

Dam 


Figure  5.  Tualatin  River  channel  bottom  elevations  as  a  function  of  River  mile.  Channel  slope 
from  RM  68  to  62  is  about  0.003  (a=0.2°)  and  the  slope  from  RM  62  to  RM  32  is  about 
0.00037(o=0.02°). 

The  Snake  River  in  Idaho  will  be  used  as  another  test  case  of  the  river  basin  model.  Figure  6  shows 
the  bottom  channel  elevation  as  a  function  of  river  mile  between  C.  J.  Strike  Reservoir  and 
Brownlee  Reservoir.  Recent  work  by  Wells  and  Berger  has  taken  the  original  W2  code  and  adapted 
it  to  this  river  system.  However,  this  implementation  required  dividing  up  a  1 10  mile  section  of  the 
Lower  Snake  River  system  into  9 1  branches  and  performing  significant  code  alterations.  Results  of 
the  new  code  development  will  be  applied  to  this  Lower  Snake  River  stretch  also  and  compared  to 
results  using  the  existing  W2  formulation. 
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Elevation,  ft  MSL 


340  360  380  400  420  440  460  480  500 

River  Mile,  Snake  River 

Figure  6.  Snake  River  channel  bottom  elevation  versus  River  mile.  Channel  slope  between  RM  450 
and  385  was  about  2.5E-4.  The  slope  between  RM  385  and  335  was  about  4.3E-4. 
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Governing  Equations  Development 


This  section  will  formally  derive  the  governing  equations  for  W2  highlighting  assumptions  and 
limitations  of  the  model  equations. 


Coordinate  System 


The  general  coordinate  system  that  will  be  used  is  shown  in  Figure  7. 


Figure  7.  Coordinate  system  for  governing  equations  (x  is  oriented  E,  y  is  oriented  N,  and  z  is 
oriented  upward). 
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Note  that  Q  is  a  vector  that  represents  the  angular  velocity  of  the  earth  spinning  on  its  axis.  The 
rotation  of  our  coordinate  system  can  result  in  significant  horizontal  accelerations  of  fluids.  This  is 
usually  restricted  to  large  water  bodies,  such  as  large  lakes  and  ocean  systems.  The  force  that 
causes  horizontal  accelerations  as  a  result  of  the  spinning  coordinate  system  is  termed  the  Coriolis 
force. 

Turbulent  Time-Averaged  Equations 

The  governing  equations  are  obtained  by  performing  a  mass  and  a  momentum  balance  of  the  fluid 
phase  about  a  control  volume.  The  resulting  equations  are  the  continuity  (or  conservation  of  fluid 
mass)  and  the  conservation  of  momentum  equations  for  a  rotating  coordinate  system  (Sabersky  et 
al.,  1989;  Cushman-Roisin,  1994;  Batchelor,  1967).  After  using  the  coordinate  system  in  Figure  7, 
applying  the  following  assumptions: 

•  incompressible  fluid 

•  centripetal  acceleration  is  a  minor  correction  to  gravity 

•  Boussinesq  approximation 

—  = - - - »  —  where  p  =  p0  +  Ap  where  po  is  a  base  value 

p  p0+Ap  pQ 

and  Ap  has  all  variations  in  p 

and  substituting  the  turbulent  time  averages  of  velocity  and  pressure  as  defined  below 

•  all  velocities  and  pressure  are  considered  the  sum  of  turbulent  time  averages  and 

1  T 

deviations  from  that  average,  i.e.,  u  =  u  +u' ,  where  u  =  —  udt  as  shown  in  Figure 

8.  The  other  terms  are  v  =  v  +  v'  ;w  =  w  +  w'  and  p  —  p  +  p'  where  the  overbar 
represents  time  averaged  and  the  prime  represents  deviation  from  the  temporal  average; 
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the  governing  equations  become  after  simplification: 


Continuity 


du  dv  dw 
dx  dy  dz 


=  0 


where  u,  v,  w  are  the  velocities  in  the  x,  y,  and  z  axes,  respectively; 


Longitudinal  momentum 


_du  _du  _du  _ 

■U  —  +  V  —  +W—-  2Qv+2  Qw  = 
dx  dy  dz  v _ i 


convective  acceleration 


Coriolis  acceleration 


viscous  stresses 


turbulent  stresses 


where  x„:  turbulent  shear  stress  acting  in  x  direction  on  the  x-face  of  control  volume  (see  Fig.  9) 
xxy:  turbulent  shear  stress  acting  in  x  direction  on  the  y-face  of  control  volume  (see  Fig.  9) 
T*z'-  turbulent  shear  stress  acting  in  x  direction  on  the  z-face  of  control  volume  (see  Fig.  9) 
Li:  dynamic  viscosity 

Q,:  component  of  Coriolis  acceleration  where 
Qz:  Q.e  sin0 

Qy:  Qe  cos0 

<{>:  latitude  of  the  earth 

dE:  rotation  rate  of  the  earth 
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Figure  9.  Sketch  of  turbulent  shear  stresses  in  x-direction. 


Lateral  momentum 

dv  _dv  _dv  _dv  _  „  _ 

—  +  u—+v  —  +  W  —  +  2Q.M  —2L2xw  = 
at  ox  ay  oz 

ydp_  drv_  £*»  dz>?  1 

p  dy  p  V  dx 2  dy2  dz 2  )  P  ^  dx  dy  dz  J 

where:  Tyx:  turbulent  shear  stress  acting  in  y  direction  on  the  x-face  of  control  volume  (see  Fig.  10) 

Xyy:  turbulent  shear  stress  acting  in  y  direction  on  the  y-face  of  control  volume  (see  Fig.  10) 

V  turbulent  shear  stress  acting  in  y  direction  on  the  z-face  of  control  volume  (see  Fig.  10) 

Qx=0 
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Figure  10.  Sketch  of  turbulent  shear  stresses  in  y-direction. 


Vertical  momentum 

dw  _dw  _dw  _ dw  _ 

+  u  -r— +  v  -r— +  w~r—  —  2Q.ru  +  20  v  =  —  e 
at  dx  dy  dz  '  J 

_1_  dp  jJ.  r d2w  d2w  d2w'  1  ^ dzzy  dz„  'j 

P  dz  p\dx2  dy2  dz2  ;  +  p  ^  dx  +  dy  +  dz  ) 

where:  za:  turbulent  shear  stress  acting  in  z  direction  on  the  x-face  of  control  volume  (see  Fig.  1 1) 

V  turbulent  shear  stress  acting  in  z  direction  on  the  y-face  of  control  volume  (see  Fig.  11) 

Xz,:  turbulent  shear  stress  acting  in  z  direction  on  the  z-face  of  control  volume  (see  Fig.  1 1) 

£2X= 0 


Figure  1 1 .  Sketch  of  turbulent  shear  stresses  in  y-direction. 

Note  that  the  turbulent  shear  stresses  are  defined  as  follows: 
t  xx  =  pu'u' 

Trv  =  pu'v'  is  the  same  as  t xx  =  pv'u ' 

rxz  =  pu'w'  is  the  same  as r„r  =  p  w'u' 
rvv  =  pW7 

rvz  -  pv'w'  is  the  same  as  x  =  p  w~'v ' 

T  ..  =  p  W  '  w  ' 

Coriolis  Effect 
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As  noted  above,  all  the  Qx  terms  are  zero  and  can  be  eliminated  from  the  y  and  z-momentum 
equations.  If  one  integrates  over  the  y-direction  (therefore  assuming  the  net  velocity  in  y  is  zero) 
and  assumes  that  the  horizontal  length  scale  is  much  greater  than  vertical  length  scale,  it  can  be 
shown  by  using  scaling  arguments  that  the  Coriolis  acceleration  forces  are  negligible  (Cushman- 
Roisin,  1994).  Hence,  prior  to  lateral  averaging,  the  Coriolis  acceleration  terms  will  be  neglected. 

Coordinate  System  Transformation 

The  coordinate  system  will  be  transformed  into  a  form  compatible  with  the  original  W2 
development  where  the  vertical  axis  is  in  the  direction  of  gravity.  In  addition,  as  shown  in  Figure 
12,  the  coordinate  system  will  be  oriented  along  an  arbitrary  slope. 


Figure  12.  General  coordinate  system  with  z-axis  compatible  with  original  derivation  of  W2  model. 
The  gravity  acceleration  is  a  body  force  that  is  then  represented  by  a  vector: 
g  =  -gVh 


where  h  is  the  surface  normal  from  the  earth’s  surface  (h  is  an  elevation  in  the  opposite  direction  to 
the  acceleration  of  gravity  vector)  and  g  is  the  acceleration  constant  (9.8  m/s  ). 

This  term  can  be  written  as  three  vector  components: 


19 


Figure  13.  Sketch  of  channel  slope  and  coordinate  system  for  W2  where  the  x-axis  is  oriented 
along  the  channel  slope. 


The  channel  slope  can  also  be  incorporated  into  the  definition  of  the  gravity  vector  if  the  x-axis  is 
chosen  parallel  to  the  channel  slope  as: 

The  channel  slope  is  defined  as  S  =  sin  ct  =  — — 

ox 


And  hence, 


dh 

gx  =~g—  =  gsma 
dh 


gz  =-g~^  =  gcosa 


The  gravity  acceleration  in  y  is  assumed  to  be  negligible  since 
channel. 


dh 

~x~  =  0  in  the  lateral  direction  of  the 
dy 
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Governing  Equations  for  General  Coordinate 
System 

After  making  the  following  simplifications: 

•  redefine  coordinate  system 

•  eliminate  Coriolis  effects 

•  neglect  viscous  shear  stresses 

The  governing  equations  become: 


Continuity 

du  dv  dw 

Longitudinal  momentum 


du  _du  _dii  _du  .  1  dp 

—  +u  —  +  v  —  +  w  —  =gsina  -—3-  ■ 
at  ox  ay  oz  v — v — '  p  ox 

' — v — '  v _  v  gravity  s _ l_v _ j 

unsteady  acceleration  convective  acceleration  pressure  gradient 


1  ( dlrr  dt  d% 


turbulent  shear  stresses 


Lateral  momentum 


dt  _dv  _dv  _dt  1  dp  \(dr,x  dxn.  dxyz' 
~^  +  Ud^  +  V  dy  +  W  dz~  pdy  + p{  dx  +  dy  +  dz 

Vertical  momentum 


dw  _dw  _  dw  dw  1  dp 

—  +  u  —  +  v  —  +  w  —  =  gcosa-  — — 
dt  dx  dy  dz  p  dz 


1  f  dx „  dx  dx..  \ 
p\  dx  dy  dz  ) 


Simplification  of  Vertical  Momentum  Equation 

Assuming  that  the  longitudinal  length  scale  is  much  greater  than  the  vertical  length  scale,  then  the 
vertical  velocities  «  horizontal  velocities.  A  result  of  this  assumption  is  that  vertical  velocities  are 
very  small  such  that  the  z-momentum  equation  becomes  the  hydrostatic  equation: 

1  dp 

— —  =  gcosa 
p  dz 
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This  assumption  prevents  the  model  from  accurately  modeling  vertical  accelerations  of  the  fluid  as 
a  result  of  convective  cooling  at  night  and  other  such  vertical  accelerations. 

Lateral  Averaging  of  3-D  equations 

The  governing  equations  above  will  be  laterally  averaged  after  decomposing  all  velocities 
and  pressure  into  a  lateral  average  and  a  deviation  from  the  lateral  average.  The  vertical 
and  longitudinal  velocities  and  pressure  are  defined  as  follows: 

u  —  u  +  u  where  u  =  -j  udy  and  B  is  the  width  of  the  control  volume 

w  =  w  +  w" 
v  =  v +v" 

p  =  p  +  p" 

The  double  overbars  represent  the  spatial  average  of  the  temporal  average  quantity.  The  double 
prime  represents  the  deviation  from  the  lateral  average  and  is  a  function  of  y.  This  is  shown  in 
Figure  14. 


Figure  14.  Lateral  average  and  deviation  from  lateral  average  components  of  longitudinal  velocity. 

These  definitions  are  substituted  into  the  turbulent  time-average  governing  equations  and  then 
laterally  averaged.  The  y-momentum  equation  is  neglected  since  the  average  lateral  velocities  are 
zero,  i.e.,  v  =  0. ,  and  cross  shear  stresses  that  contribute  to  vertical  mixing  will  be  computed  from 


the  analysis  of  wind  stress.  The  equations  that  remain  are  the  continuity,  x-momentum,  and  z- 
momentum  equations. 

Continuity 

The  continuity  equation  becomes  after  substituting  the  above  velocity  components  and  laterally 
averaging 


d(u+u ")  d(v  +  v")  d(w  +  vQ 

dx  dy  dz 

The  lateral  average  of  a  double  primed  variable  is  by  definition  zero,  i.e., 

=  1  -v2 

u"  =  —  \u"dy  =  0 

B  vl 


Also,  note  that: 


d(v  4-  v")  1  yfd(v  +v")  _  (v  +v") 

~~B  J  “  y~ 


dy 


dy 


B 


B 


=  <? 


where  q  is  defined  as  the  net  lateral  inflow  per  unit  volume  of  cell  [T1] 


d(H  +  u")  _  l_yrd(u  +  u")  7  _  1  fdtf  ^  ,  J_f  du 

B  J 


dx 


dx 


i  rdu  1  row  1  d  'r_  1 

^  =  B  J  A  +  «  J  ^  =II  1  ' = 

vl 


B  \  dx 

Vl 


B  dx 


and 


d(W  +  w")  1  -vrd(w  +  w")  _  1  ydw  1  ydw"  _  1  d  r  = .  1 


A 


-if 


dl 


\  1  f  <7VV  I  f  <7W  1  U  f  — 


Combining  terms,  the  continuity  equation  becomes 
dBu  dBw 

-&+ir-qB 

Longitudinal  momentum 


The  laterally  averaged  x-momentum  equation  is  more  easily  simplified  by  writing  it  in  conservation 
form  (this  can  be  verified  by  using  the  continuity  equation  with  the  x-momentum  equation). 
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d(u+u")  d(u+u")(u  +  u")  d(v  +  v")(u  +  u")  d(w  +  w")(u  +  u") 

dt  dx  +  dy  +  * 

.  1  d(p  +  p")  lfdzxx  dt  dz  \ 

gsma-— - -f— + - -f +  — -  +  — ^ 

p  ox  P\  OX  ay  dz  J 


Each  term  in  this  equation  can  be  simplified  as  follows  (note  that  the  spatial  average  of  any  double 
primed  variable  goes  to  zero  by  definition): 


Unsteady  acceleration  term 


d{u+u")  1  p(u+u")  J  1  }fdu  ,  1  Y  du"  ,  1  d  Y_ .  1  d  Y  1  dRlf 

— — b>~ 'sr~dy=z)  k<*>— 

vl 


_  i  rau  1  [du  _  1  d  r __ ,  1  d  r  1 

“«J  J“  dy=l 


B  dt 


Convective  acceleration  term 


d(u+u")(u+u")  1  'f d(u  +u"){u+u")  ,  1  Y duu,  1  y2 duu"  1  yrdu"u" 

*  =b[  * - dy=i{-^rdy+B[-^rdy+^i—dy 


juudy 


1  d  y 

u"u"dy  = 


1  dBuu  1  d  Y 

^~&~+i&{u"u''dy 


Bit  & 


dispersion  term 


Similarly  for  the  other  two  terms: 

d(u  +u"){w  +  w")  1  dBuw  1  d  'f 

a  =  t ~  *  1  a  rw"dy 


dispersion  term 


d(u  +  u")(v  +  v") 


—  ,  ''  ,//  //  //|  o 

—  MV  -LI  V  |v|  =0 


Gravity  term 


gsina  =  — Jg  sin  atfy  =  -j-(g  sin  a)  Jt/y  =  gsina 


24 


Pressure  gradient  term 


d(p  +  p")  lyfd(p  +  p")  ,  If }dp  ,  if  dp"  ,  1  <9  f =J  ,  1  <9  f  //; 

-*  ■  i { — *— * -  *  J  *  *+ « hr*  -  i** + **  J '  *  “ 
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or  the  above  equation  can  be  written,  assuming  that  the  derivative  of  the  lateral  average  pressure 
gradient  in  the  x-direction  is  not  a  function  of  y: 


+  p") 
dx 


l  p(p  +  /U) 

5v.  * 


i  dp  f  r  l  f  <9p' 

*-:«*J*+»/- 


<9c 


■Jy 


B  dx 


v2 


1  <9  Y 


£ 

<2r 


Shear  stress  term 


dxxx  dr„  dx 

- ±L  _| - -  ^ - — 

dx  dy  dz 

1  (9 


+  -J~By  +  -  f -~-dy  =  -4~j  rjy  + 

Bl  dx  y  Si  dy  y  B  *,  dz  y  Rd* J  “ 


l  d  f  ,  i  a  r  i 


v 


<9bt„  dBx xx  dBrx, , 

- 1 - 1 - 

dx  dy  dz 


B  dx 

^  '  ( dBxrr  dBx  ^ 

+ 


B 


dx  dz 


Then  collecting  all  terms  and  neglecting  all  dispersion  terms,  the  final  x-momentum  equation  is 
then  after  simplification: 


dBu  dBu  u  dBu  w 

dt  dx  dz 


Bg  sin  a  - 


Bdp_ 
p  dx 


+ 


*0*2 

dz 


Summary  of  Laterally  Averaged  Equations 

In  the  development  of  W2  in  Cole  and  Buchak  (1995),  the  lateral  average  terms  were  represented 
by  uppercase  characters,  such  that  ui  =  U ,  w  =  W ,  and  p  =  P .  The  shear  stress  terms  will  be 
assumed  to  be  lateral  averages  and  the  double  overbars  will  be  dropped  for  convenience.  Making 
these  simplifications,  the  governing  equations  become 

Continuity 


BUB  bwb 

dx  dz 


qB 


Longitudinal  momentum 


auB  auuB  awuB  #ap  ibbt 

at  ox  dz  p  dx  p  dx 


Vertical  momentum 

1  dP 

— — =  gcosa 
p  dz 


Now  we  have  three  equations  and  three  unknowns:  U,  W,  and  P. 


1  aBrx, 

p  dz 
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Simplification  of  the  Pressure  Term 


The  z-momentum  equation  reduces  to  P  =  Pa  +  g  co saj'pdz  after  integration  from  a  depth  z  to 

the  water  surface  defined  as  z=r\.  Pa  is  the  atmospheric  pressure  at  the  water  surface  (see  Figure 
15). 


Zl  ^surface  "H 


Figure  15.  Illustration  of  layout  for  simplification  of  pressure  term. 


This  equation  for  pressure  is  now  substituted  into  the  x-momentum  equation  and  simplified  using 
Leibnitz  rule.  The  pressure  gradient  term  in  the  x-momentum  equation  then  becomes: 


J dP 
p  dx 


1  dP 

ci 

p  dx 


+  g 


dp 


cos  cc—- 
dx 


gcosa 

P 


The  first  term  on  the  RHS  is  the  atmospheric  pressure  term  that  accounts  for  accelerations  due  to 
atmospheric  pressure  changes  over  the  water  surface.  The  second  is  the  barotropic  pressure  term 
that  accounts  for  accelerations  due  to  water  surface  slopes,  and  the  third  is  the  baroclinic  pressure 
term  that  accounts  for  accelerations  due  to  density  differences. 

In  W2,  the  atmospheric  pressure  term  is  assumed  to  be  zero  and  is  neglected.  This  implies  that  for 
long  systems  during  severe  storms  the  model  will  not  be  able  to  account  for  accelerations  on 
account  of  atmospheric  changes.  For  a  large  physical  domain,  variations  in  meteorological  forcing 
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may  be  significant.  This  is  discussed  in  Variability  in  Meteorological  Forcing.  The  pressure  term 
then  becomes  with  this  simplification: 


1  dP  dri 

---=gcosa-- 


geos  a 
P 


The  revised  form  of  the  x-momentum  equation  is  then: 


5UB  auUB  dWUB 

+  T -  +  ” - 

8t  dx  dz 


n  .  ndr\  gcosaB  rdp 

gBsma  +  gcosocB—-- -  -$-dz  + 

ax  p  ;  dx 


\  d#7,, 

p  dx 


1  3Btx, 
p  dz 


Effectively,  pressure  has  been  removed  from  the  unknowns  by  combining  the  z-momentum  and  x- 
momentum  equations,  but  we  have  added  r|  as  an  unknown. 


Free  Water  Surface  Equation 

This  equation  is  a  simplification  of  the  continuity  equation.  The  continuity  equation  integrated  over 
the  depth  from  the  water  surface  to  the  bottom  is  called  the  free  water  surface  equation.  Figure  15 
and  Figure  16  are  definition  sketches  for  the  W2  cell  layout  without  and  with  a  channef  slope, 
respectively. 


CE-QUAL-W2  coordinate  system:  a=0 


g 


28 


CE-QUAL-W2  coordinate  system:  a>0 


Figure  17.  W2  coordinate  system  with  finite  channel  slope. 
The  continuity  equation  is  integrated  over  the  depth  as  follows: 


fdUB,  (9WB,  }  nJ 

i^dz+  i~^rdz  =  )qBdz 


dx  -  ■  i  dz 
n  n 


The  first  term  can  be  expanded  as  follows  using  Leibnitz’s  rule: 


r3UB  a  r  an  ,  on _ , 


5  f 


dh 


dr] 


dx 


dx 


dx 


dx 


The  integral  of  the  vertical  flow  rate  over  z  relates  to  changes  in  water  surface  elevation  as  shown 
below: 


r  awB  .  . 

=WB\h-WB\n 


dh  dh 

where  + 


Combining  these  terms  together,  the  free  surface  equation  becomes: 


dx 


dr] 

dx 


dh 


dh 


dr] 


dr\ 


n 

=  J  qBdz 


Canceling  out  terms  and  applying  the  no-slip  boundary  condition  that  Uh  is  zero: 

ljuBdz-B,f  =  ]qBdz 


or 

*1  n 

where  Bn  is  the  width  at  the  surface. 

Equation  of  State 

The  density  must  be  known  for  solution  of  the  momentum  equations.  The  equation  of  state  is  an 
equation  that  relates  density  to  temperature  and  concentration  of  dissolved  substances  and  is  given 
as  follows: 

P  =  f(Tw ,  ‘J’tds  ,  Oss) 

where  f(Tw,0TDs,  Oss)=density  function  dependent  upon  temperature,  total  dissolved  solids  or 
salinity,  and  suspended  solids. 

Hence,  the  temperature,  total  dissolved  solids,  and  suspended  solids  must  be  known  and  are 
determined  from  the  water  quality  model. 

Summary  of  Governing  Equations 

Table  1  shows  the  governing  equations  after  lateral  averaging  for  a  channel  slope  of  zero  (original 
model  formulation)  and  for  an  arbitrary  channel  slope.  Parameters  used  in  Table  1  are  illustrated  in 
Figure  17. 


Table  1.  Comparison  of  governing  equations  for  W2  with  and  without  channel  slope. 


Equation 

Existing  governing  equation  assuming  no  channel 
slope 

Governing  equation  assuming  an  arbitrary  channel  slope 

X- 

momentum 

3ub  3uub  3wub 

9t  9x  dz 

dT)  gB\  dp 

8 

1  1  SBtxj 

- + - 

p  dx  p  dz 

9  UB  9UUB  9  WUB 

1-  4-  — 

d  t  dx  dz 

dr\  g  cos  aB  }  dp 

gB  sin  a  +  g  cos  aB  -  I  ^  dz  + 

dx  p  ^  dx 

1  dBr 1  3Brx- 

- + - - 

p  dx  p  dz 

z-momentum 

1  dP 
°~8~  pdz 

1  dP 

0=  geos  a  - — — 
p  dz 

free  surface 
equation 

^\uBdz-\qBdz 

n  ti 

Note:  U,W:  horizontal  and  vertical  velocity 
P:  pressure 

tx,tz:  lateral  average  shear  stress  in  x  and  z 
rp  water  surface 


B:  channel  width 
g:  acceleration  due  to  gravity 
p:  density 


a:  channel  angle 


grav  ity 


channel  slope  =  S  0  =  sin  a 


W~M~D  Datum 

Figure  18.  Definition  sketch  for  channel  slope  (exaggerated  slope). 


Determination  of  Turbulent  Shear  Stresses 

In  order  to  solve  the  equations  above,  the  following  terms  must  also  be  evaluated: 

•  shear  stress  at  water  surface  from  wind 

•  shear  stress  at  boundary  (boundary  shear) 

•  algorithm  for  transmitting  shear  stress  vertically, 

•  algorithm  for  transmitting  shear  stress  longitudinally,  xxx 

Surface  shear  stress 

The  shear  stress  at  the  water  surface  is  defined  as: 

^  =  cDp„(w'-«,)2  =  c0p„(iv,)2 

where  xs:  surface  shear  stress  at  water  surface 
us:  surface  velocity  in  water 

Wh:  wind  velocity  measured  at  a  distance  h  above  water  surface  in  direction  of  shear 
CD:  drag  coefficient 
pa:  air  density 
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Note  that  this  relationship  leads  to  the  “3%  rule”  for  surface  currents: 


~  CD  Pa  (Wh  wt)  —  CDpwus  i  if  CD  ~  CD 


then  us  ~  0.03Wh 

' - » - ' 

3%  rule 


Usually  the  drag  coefficient  is  a  function  of  the  measurement  height,  h,  above  the  water  surface. 
Most  drag  coefficient  formulae  have  been  determined  based  on  a  10  m  wind  speed  measurement 
height.  If  wind  speeds  are  taken  at  other  measurement  heights,  for  the  shear  stress  calculation,  these 
should  be  corrected  to  10  m. 


The  windspeed  is  a  function  of  measurement  height.  To  correct  the  measurement  height  to  an 
elevation  z,  the  following  approach  is  used. 


Assuming  a  logarithmic  boundary  layer: 


Wt 

Wr. 


inf—; 

Zo 

Inf—; 

Zo 


where  Wz:  desired  wind  speed  at  elevation  z 
Wz]:  known  wind  speed  at  height  Z\ 

Zo:  wind  roughness  height  (assume  0.003  ft  for  wind  <  5  mph  and  0.015  for  wind  >  5 
mph,  range  0.0005  to  0.03  ft) 


This  term  can  then  be  used  to  compute  the  surface  stress  in  the  direction  of  the  x-axis  and  the  cross¬ 
shear  (the  cross-shear  term  will  be  used  in  the  turbulent  shear  stress  algorithm)  as  follows: 


cos(0,-02) 
s  Q  p,  iv;!  sin<0  -0,) 


where  t**:  surface  shear  stress  along  x-axis  due  to  wind 

xwy:  surface  shear  stress  along  lateral  direction  due  to  wind 
0i :  wind  orientation  relative  to  North,  radians 
02:  segment  orientation  relative  to  North,  radians 
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The  drag  coefficient,  Cd,  is  defined  in  W2  as: 

For  Wh<l  m/s,  CD  =  0.0 

For  1<  Wh<15  m/s,  CD  =  0.0005(Wh)°'5 

For  Wh>  15  m/s,  CD  =  0.0026 

These  formulae  are  based  on  a  10  m  measurement  height.  Note  that  for  a  river  basin  model,  the 
wind  speed  may  vary  significantly  over  the  study  area.  This  issue  is  discussed  in  the  section  on 
Variability  in  Meteorological  Forcing. 


Bottom  shear  stress 

The  shear  stress  is  defined  along  the  bottom  of  each  cell  or  for  each  cell  in  contact  with  side  walls 

or  channel  bottom  as: 

rb=^fup\ 

c 

where  C  is  the  Chezy  friction  coefficient 
U  is  the  longitudinal  velocity 
pw  is  the  density  of  water 

The  Chezy  coefficient  is  related  to  the  Manning’s  friction  factor  as 

C  (for  SI  units  only)=  (l/n)R1/6 

where  n:  Manning’s  friction  factor 
R:  hydraulic  radius 


The  shear  stress  for  bottom  friction  will  be  more  important  for  sloping  river  channels  since  the 
primary  balance  is  between  gravity  and  bottom  friction.  The  bottom  friction  coefficient  will  be 
variable  segment  by  segment. 

Vertical  shear  stress 

Vertical  shear  stress  in  W2  is  defined  as: 

dU  dU 

p  ~  V'“rh“lm'  dz  ~  dz 

This  shear  stress  term  also  includes  the  contribution  to  the  shear  stress  from  surface  waves  induced 
by  the  wind.  The  wind  can  produce  waves  that  produce  decaying  motions  with  depth  (Lighthill, 
1978)  as  shown  in  Figure  22. 
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wind  stress 


Figure  22.  Decay  of  wind  shear  with  depth. 


This  addition  has  been  critiqued  by  some  (Chapman,  1989)  as  being  unnecessary  since  the 
expressions  for  the  wind  shear  stress  implicitly  account  for  the  effect  of  wave  induced  stress  by  the 
wind  shear. 


Regardless,  the  total  longitudinal  shear  stress  for  a  layer  is  defined  in  W2  as  having  contributions 
from  interfacial  velocity  shear,  wind  wave  generated  shear,  and  friction  shear  along  boundaries: 


£s L 

P 


9z  p 


-2  kz 


P 


where  Twx  is  the  longitudinal  wind  shear  at  the  surface  (see  above) 

4x2 

k=  wave  number  = - =- 

gT; 

Tw=  wind  wave  period  (empirical)  =  6.95 E  —  2  F0'233  |W,|°'5j’4 

F=  fetch  length,  m. 

Vertical  eddy  viscosity 


The  turbulent  eddy  viscosity  was  conceptualized  by  Prandtl  as: 


Y  turbulent 


=  c- 


dU_ 

dz 
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where  t  is  defined  as  the  mixing  length  and  can  be  interpreted  as  being  proportional  to  the  average 
size  of  large  eddies  or  the  length  scale  of  a  turbulent  eddy.  This  length  is  a  function  of  distance 
from  a  boundary  or  wall  since  the  eddy  sizes  vary  as  a  function  of  distance  from  a  boundary.  The 
goal  in  most  turbulence  models  is  the  determination  of  the  mixing  length  as  a  function  of  position  in 
the  fluid. 

In  the  formulation  in  W2,  the  mechanism  for  transporting  the  wind  stress  at  the  surface  is  based  on: 


Az=  vertical  eddy  viscosity  =  K— 


(dU\-  (dV\ 
+ 


dz 


dz 


1/2 


-CRi 


where: 


Ri:  Richardson  number  = 


dp 
8  dz 


f— T 

dz 


k  is  the  von  Karman  constant  =  0.4 

C  is  an  empirical  constant  taken  as  1 .5 

1,  a  vertical  length  scale,  is  chosen  as  vertical  cell  thickness. 


Hence,  this  formulation  is  a  mixing  length  formulation  that  is  decreased  or  increased  based  on  the 
Richardson  number.  The  Richardson  number  accounts  for  the  impact  of  density  stratification  on 
transfer  of  momentum  between  fluid  parcels.  In  regions  where  there  is  no  stratification,  Ri=0,  and 

dp 

the  exponential  term  is  1 .  For  regions  where  there  is  strong  stratification  (or  as  —  — »  00 ),  the 
Richardson  number  becomes  large  and  the  exponential  term  approaches  0. 

The  term  in  the  above  formulation  uses  the  lateral  velocity  because  even  winds  blowing  at  right 
angles  to  the  model  cell  may  cause  vertical  transfer  of  momentum,  or  mixing  between  fluid  layers. 
Hence,  this  term  allows  for  increasing  the  transfer  of  stress  vertically  in  the  fluid  as  a  result  of 
lateral  currents. 

In  the  longitudinal-vertical  model,  the  lateral  velocity,  V,  and  its  gradient,  3V/dz,  are  due  to  the 
lateral  component  of  wind  wave  motion  and  are  assumed  to  be  zero  when  averaged  laterally,  but  not 
necessarily  the  term  (3V/3z )’.  It  is  assumed  that  cross  wind  shear  xwy  generates  lateral  wave 
components  and  decays  exponentially  with  depth,  z,  such  that 

Tyz=xwy  exp(-2kz) 

where  xwy  is  the  lateral  wind  shear  at  the  surface  defined  above. 
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Then  using 

—  =  A  — 

p  ~  z  dz 

the  lateral  velocity  gradient  squared  becomes 


rdV' 

2 

T,v  exp(-2 kz) 

k  dz  J 

PAz 

The  final  equation  for  the  vertical  eddy  viscosity  is  then 


A,  =  k\  — 


fx  e-2kzV 

t  Wy 


P  A 


„(-CRi) 


Z  J 


The  above  equation  is  implicit.  In  the  model,  this  equation  is  explicit  since  the  value  of  Az  in  the 
lateral  wind  shear  term  is  used  from  the  previous  time  step.  Az  is  never  less  than  the  molecular 
kinematic  viscosity  for  water.  This  formulation  is  subject  to  further  research  and  discussion  and  is 
vital  in  establishing  flow  patterns  within  the  domain. 

In  the  river  basin  sections,  the  above  formulation  for  vertical  eddy  viscosity  may  or  may  not  be 
adequate.  Additional  formulations  for  the  vertical  eddy  viscosity  for  both  vertically  well-mixed 
systems  and  continuously  stratified  systems  will  be  tested  as  shown  in  Tables  2  and  3.  Since  these 
formulations  vary  significantly,  the  theory  and  literature  values  will  be  used  as  a  guide  to  determine 
the  appropriate  function  after  model  testing  (this  follows  the  recommendation  of  Shanahan  and 
Harleman,  1982).  Numerical  results  of  the  velocity  field  will  be  compared  to  analytical  solutions 
for  2-D  open  channel  flow  at  steady-state  with  simple  turbulence  closure  hypotheses.  The  code  will 
be  tested  for  stability  and  robustness  of  the  solution  scheme. 


Table  2.  Typical  vertical  eddy  viscosity  (Av)  formulations  for  shallow  systems  that  are  vertically 


well-mixed  where  x  —  A  — . 

v'  dz 


Equation 

number 

Equation 

Comments 

Reference 

1 

ii 

\h-z-  0.1) 

h  J 

3/4 

Avo  is  the  eddy  viscosity  at  the 
surface,  zero  at  bottom, 
maximum  at  surface 

Fjeldstad  in  Neumann 
and  Pierson  ( !  966) 

2 

Av  =  KUzZ 

( 

1-1 

\  h  y 

parabolic  shape  (zero  at  top 
and  bottom) 

Engel  und  (1978) 

3 

Av  =KU.(1-Z)  7 
\hj 

parabolic  shape  (zero  at  top 
and  bottom) 

Koutitas  (1978) 

4 

i 

s 

ii 

< 

h  j 

1/2 

+  4, 

minimum  value  at  bottom, 
maximum  at  water  surface 

Liggett  ( 1 970) 
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Equation 

number 

Equation 

Comments 

Reference 

5 

Av  =  KUtZ  for  near  —  surface 

Av  -  Kuih  (h  -  z )  for  near  -  bottom 

u*b  is  the  bonom  shear 
velocity,  zero  at  surface  and 
bottom 

Madsen  (1977) 

6 

minimum  at  bottom  and 
maximum  at  surface 

Witten  and  Thomas 
(1976) 

7 

Av  =  —  ~J 

maximum  at  surface  and  zero 
at  bottom,  Av0  is  Ku*h 

Thomas  (1975) 

8 

(  z  V 

A-  =  A~V~wh) 

Maximum  at  surface  and 
minimum  at  z=h,  AVo=C]\vh, 
W—  l-(Avh^Avo)  »  Cj  is  a 
constant 

Lindijer  (1979) 

Note:  z  is  zero  at  the  surface  and  is  positive  downwards,  hence  at  the  bottom  z=h,  where  h  is  the 
total  depth,  u.  is  the  shear  velocity,  K  is  von  Karman  constant.  AVh  is  the  bottom  value  of  eddy 
viscosity. 


Table  3.  Typical  vertical  eddy  viscosity  (Av)  formulations  for  systems  that  are  continuously 

du 

stratified  where  Tz  =  Av  — . 


Equation 

number 

Equation 

Comments 

Reference 

1 

A,  =  Ajl  +  PRi)“ 

Avn  is  the  eddy  viscosity  at  the  surface, 
Ri  is  the  Richardson  number, 

dp 

r.  8  dz 

Kl  —  /  ?  ,  a  and  p  were 

(  du) 

empirical  coefficients* 

Munk  and  Anderson  (1948) 

2 

(  A  Y' 

a^Wk) 

Avo  is  the  eddy  viscosity  at  the  surface, 

D  ghAp 

Ro  is  defined  as  K  —  ?  ,  ai 

pul 

and  pi  were  empirical  coefficients* 

French  ( 1 979) 

3 

Av  =  Am  exp (~p2Ri) 

Avo  is  the  eddy  viscosity  at  the  surface, 

p2  is  an  empirical  coefficient* 

Mamayev  as  quoted  in  French 
(1985) 

The  range  of  coefficients  values  are  shown  from  French  (1985)  to  be:  for  a=-l ,  (3  ranges  from  2.5  to 
30.3;  a=-2,  P  =9.8;  a=l,  (3  =-3.3;  a=-0.5,  (3  ranges  from  10  to  30;  a, =0.747,  pi=0.31;  (Xi=0.379, 
Pi  =0.062;  and  for  model  3,  p2=0.4; 


Longitudinal  turbulent  shear  stress 

The  longitudinal  turbulent  shear  stress  is  defined  as: 

_  dU  dU 
~p  ~  V'"rhl'u""  rjx  ~  A*  dx 
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where  Ax  — v^^ient  and  is  the  longitudinal  turbulent  viscosity  or  the  longitudinal  eddy  viscosity.  Ax 
is  a  user-defined  constant  in  the  model. 

This  turbulence  closure  approximation  is  termed  a  zero-order  closure  model  since  no  further 
equations  are  necessary  to  solve  for  the  transmission  of  longitudinal  shear  stress  within  the  fluid. 
This  term  is  usually  of  very  low  magnitude  except  in  areas  near  boundaries,  like  at  the  face  of  a 
dam,  where  the  longitudinal  velocity  goes  to  zero. 

Validation  of  Velocity  Field 

The  velocity  predicted  by  W2  in  river  channels  with  sloping  bottoms  will  be  compared  to  the 
following  theoretical  approaches  for  open-channel  flow: 

•  steady-state,  depth  averaged  velocity  will  be  compared  to  Chezy  or  Manning’s  Equation,  i.e., 

1  h  1 

— - -  J  Udz  =  C^fRS  =-RmS'n 

( h-V)J„  n 

where  C  is  the  Chezy  coefficient 

n  is  the  Manning’s  friction  factor 
S  is  the  channel  longitudinal  slope 
R  is  the  hydraulic  radius 

The  model  will  be  run  at  steady-state  and  the  predictions  of  velocity  will  be  compared  to  the  Chezy 
and  Manning’s  Equation 


•  assuming  a  turbulent  vertical  eddy  viscosity  of  Av  =  kU,z 


1  — - 


h-ri 


,  the  velocity  profile 


varies  vertically  as  U(z)  =  ■ 


U. 


K 


1  +  ln- 


h-rj 


where  U*  is  the  shear  velocity 

k  is  the  von  Karman  constant 

This  also  is  a  steady-state  model  of  the  vertical  variation  of  velocity  in  a  sloping  channel.  The  W2 
model  results  will  be  compared  with  this  formulation,  as  well  as  others  derived  for  different 
assumed  functional  distributions  of  Av. 


Variability  in  Meteorological  Forcing 

The  river  basin  model  will  be  extended  over  large  spatial  scales  that  will  probably  include 
significant  meteorological  data  variability.  The  model  currently  requires  the  following 
meteorological  data: 
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•  air  temperature 

•  dew  point  temperature  (or  relative  humidity) 

•  wind  speed  and  direction 

•  cloud  cover  (or  measured  solar  radiation  and  a  back-computed  cloud  cover) 

Because  of  the  implications  of  spatially  variable  meteorological  data,  the  following  aspects  will  be 
evaluated: 

•  allowance  for  different  meteorological  input  files  for  different  model  branches 

•  in  river  basins  and  in  some  reservoir  systems,  shading  by  canyons  or  vegetation  in  parts 
of  upper  river  basins  may  also  be  important.  The  inclusion  of  a  simple  segment-by¬ 
segment  shading  function  that  only  adjusts  short-wave  solar  radiation  should  be 
included.  If  a  canopy  formed  over  the  segment,  then  adjustments  would  also  have  to  be 
made  on  atmospheric  radiation  inputs.  As  a  first  approximation,  this  latter  case  will  be 
neglected. 

•  if  one  meteorological  file  is  provided  with  the  elevation  of  the  meteorological  station, 
the  air  and  dew-point  temperature  and  short-wave  solar  radiation  would  be  adjusted 
according  to  elevation  of  the  system  within  the  model  and  then  allow  the  user  to  make 
branch-by-branch  adjustments  to  wind  sheltering  and  cloud  cover 

Air  temperature 

The  adjustments  to  air  temperature  would  follow  the  approach  of  Singh  (1992),  where  the  average 
temperature  (dry  bulb)  decreases  with  altitude  (assuming  an  adiabatic  process)  at  0.7°C  per  100  m. 
For  example,  in  a  modeling  study  of  Wahiawa  Reservoir  (Wells  and  Berger,  1997),  the  reservoir 
was  located  at  an  elevation  of  256.6  m  NGVD  and  the  airport  meteorological  data  were  at  an 
elevation  of  2.1  m  NGVD.  The  air  temperature  difference  between  the  stations  would  then  be  about 
1 ,8°C.  A  correlation  study  of  air  temperatures  was  performed  between  Honolulu  and  Kunia  and 
verified  that  this  temperature  correction  was  correct  within  0.1  °C. 


Dew-point  temperature 


The  air  temperature,  though,  is  not  the  only  meteorological  data  parameter  that  is  affected  by 
changes  in  altitude.  The  dew  point  temperature  also  is  a  function  of  the  air  temperature  and  relative 
humidity.  The  following  correlation  between  air  temperature,  dew  point  temperature,  and  relative 
humidity  from  Singh  (1992)  was  used  to  estimate  the  reduction  in  dew  point  temperature  between 
different  locations: 

j  _  TJew  s  (1455  +  0.1 14D(1  -  RH)  +  (2.5  +  0.0070(1  -  RH )3  +  (15.9  +  0.1 177X1  -  RH)'* 

where  T:  air  temperature  in  deg  C 

Tdew:  dew  point  temperature  in  deg  C 
RH:  relative  humidity  as  a  fraction 

With  a  known  relative  humidity  and  a  known  air  temperature  at  the  new  elevation,  the  dew  point 
temperature  at  the  new  elevation  can  be  determined.  This  equation  generally  leads  to  a  dew-point 
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temperature  reduction  that  is  of  the  same  magnitude  as  the  dry  bulb  temperature,  but  somewhat  less 
based  on  relative  humidity. 


Alternative  solar  radiation  calculation 

The  current  short-wave  solar  radiation  algorithm  in  W2  does  not  correct  for  altitude  of  the  water 
body  nor  for  atmospheric  conditions.  A  different  algorithm  can  be  included  as  a  model  option  that 
incorporates  short-wave  solar  radiation  changes  with  elevation,  as  well  as  with  atmospheric 
turbidity.  This  algorithm  has  been  used  by  Brown  and  Barnwell  (1987)  in  QUAL2E  and  the 
Environmental  Laboratory  (1986)  in  CE-QUAL-R1.  This  algorithm  includes  an  additional  tuning 
factor  called  the  atmospheric  attenuation  factor,  TURB.  This  factor  accounts  for  scattering  of  short¬ 
wave  solar  radiation  by  atmospheric  dust  and  moisture. 

To  illustrate  the  differences  in  the  two  algorithms.  Figure  23  shows  a  comparison  of  the  current  W2 
solar  radiation  formulation  and  the  algorithm  used  in  QUAL2E  and  CE-QUAL-R1  at  various 
elevations  for  an  atmospheric  attenuation  factor  of  0.25  and  a  dew-point  temperature  of  15°C. 


30-  Apr  29-  Jjn  28-  Aug  27-Ocf  26-Dec 


Figure  23.  Comparison  of  solar  radiation  formulae  for  daily  average  clear-sky  conditions  using 
current  W2  model  equation  and  those  used  in  QUAL2E  and  CE-QUAL-R  L 
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Variability  in  Reaeration  Formulae 

Since  the  river  basin  model  will  encompass  areas  that  are  dependent  on  boundary  shear  or  on  wind 
stress  for  turbulence,  the  reaeration  formula  for  these  systems  must  be  variable.  In  the  following 
two  sections,  formulae  for  reaeration  as  a  function  of  wind  speed  and  as  a  function  of  boundary 
shear  are  presented. 

Theoretical  reaeration  formulae  as  function  of  wind  speed 

Theoretical  reaeration  formulae  for  wind  speed  formulations  in  units  of  day'1  for  Ka  and  m/day  for 
K^are  shown  in  Table  4. 


Table  4.  Selected  wind  speed  formulations  for  Ka  and  KLat2o°c.- 


f  Equation 

Comments 

Reference 

K,  0.8 64 W 

K  =  — = - 

“  H  H 

W  is  in  m/s  measured  at  10  m  above 
water  surface,  H  depth  in  m 

Broecker  et  al 
(1978) 

K,  aWp 

i  JS  L 

"  H  H 

oc=0.2.  (3=1.0  for  W<3.5  m/s;  a=0.057, 
(3=2.0  for  W>3.5  m/s;  where  W  is  a 
daily  averaged  wind  speed 

Geldaetal  (1996) 

K,  0.728W0  5  -  0.3 1 7W  +  0.0372 W2 

Banks  and 

Herrera  ( 1 977) 

L“  ~  H  "  H 

! 

i 

I 

K.  0.0986W164 

a:  =  — = - 

a  H  H 

Wanninkhof  et  al. 
(1991) 

K,  0.728W0  5  -  0.3 1 1W  +  0.0372 W2 

ts  L 

U  is  the  mean  estuary  tidal  velocity  in 
m/s.  this  formula  combines  the  effect  of 
wind  from  Banks  and  Herrera  (1977) 
and  estuary  tidal  flow 

Thomann  and 

Fitzpatrick  (1982) 

~  H  ~  H 

Vf7  I 

+  3.93-^t  | 

i  i 

Do  2 

Kl  (200-  60W05)10-6 

K“=~h~  h 

Dc>2  is  the  molecular  diffusivity  of 
oxygen  (at  20°C  is  2. 1 E-9  m2/$ 

Chen 

i . — - - - - - - - 

Kl  0.362VW  CC  . 

K  =  — —  = -  W  <  5  5m  /  s 

a  H  H 

Kl  0.0277W2  . 

a  H  H 

Banks  (1975) 

K,  a+bW 

H 

W  in  m/day,  a  from  0.005  to  0.01  m/day,  \  Baca  and  Arnett  j 
b  from  1 E-6  to  1 E-5  m'1  j  (1976) 

\ 

rs  s  a  A  t  oon 'l  recommended  form  for  WQRSS  model  j  Smith  (1978) 

K,  0.64  +  0.1 28VV  !  1 

K  =  —  = -  1  1 

“  H  H  i  i 

43 


{Equation 

Comments 

Reference  j 

K,  0.156VF063 

Ka  =  =  W  <  Aim  /  5 

ti  ti 

K,  0.0269 W19 

K,=-rr  = - “ -  W  >  4.1m/  5 

ti  ti 

Liss  (1973) 

„  Kl  0.0276W2 

k°=h=  h 

Downing  and 

Truesdale  (1955) 

\  „  Kl  0.0432W2 

K-=ir=  H 

Kan  wisher (1963) 

„  Kl  0.3  \9W 
k-=h=  h 

Yu,  Tuffy,  and 
Lee (I 977) 

„  Kl  0.398 

Ka=  —  =  W<l.6m/s 

H  H 

„  Kl  0.1 5  W2 

!  Ka  =~~  = - - —  W>l.6m/s 

ti  ti 

Weiler(I974) 

„  Kl  0-5  +  0.05W2  1 

k-=h=  h 

as  used  in  current  W2  model 

Cole  and  Buchak 
(1995) 

Figure  24  shows  how  some  of  these  formulations  vary  with  wind  speed  in  the  range  of  the  winter 
wind  velocities  for  Portland. 
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Variation  of  KL  with  wind  speed  at  20C 


^3" 


ep/iu  ‘nx 


Figure  24.  Variation  of  KL  with  wind  speed  at  20  °C 


Theoretical  formulae  as  a  function  of  boundary  shear 


Table  5  shows  some  of  the  more  common  reaeration  formulae  as  a  function  of  boundaiy  shear. 


Table  5.  Reaeration  formulae  as  a  function  of  boundary  shear  at  20°C. 


Equation 

number 

Equation 

Comments 

Reference 

1 

/  \  1/2 
,  Kl  (do,v) 

a  R  HV1 

D02  is  tlie  molecular  diffusion 
coefficient  for  water  at  20oC  it  is 
8.1  E"5  fP/hr,  U  and  H  are  average 
velocity  and  deptli  of  channel; 
D02-l  .9 1  E3( 1 . 03 7)r'20  where  D  is 
in  ft2/day 

O’Connor  and  Dobbins 
(1958) 

2 

to  § 

II 

II 

U  in  fps  and  H  in  ft  and  Ka  m  day'1 

Churchill,  Elmore  and 
Buckingham  (1962) 

3 

Ka  =  0.88 US  for  10  <  0  <  300 cfs 
Ka  =  1 MJS  for  1  <Q<  10  cfs 

S  is  slope  in  ftAmle,  U  is  velocity  in 
fps,  Ka  is  in  day'1 

Tsivoglou  and  Wallace 
(1972) 

4 

„  Kl  2L6U 067 
a~  H  ~  HlS5 

U  in  fps,  H  in  ft 

Owens  et  al.  (1964) 

5 

K  =i=25“'(1+F»S) 
a  H  H  y  ’ 

u*  is  tlie  shear  velocity  [(HSg)0  5],  S 
is  tlie  slope  of  energy  grade  line,  F 
is  Froude  number  [u*/(gH)° 5] 

Thackston  and  Krenkel 
(1966) 

6 

„  Kl  7.62U 
a~  H  ~  Hv 33 

U  is  in  fps  and  H  is  in  ft 

Langbien  and  Dunun  (1967) 

Temperature  dependence  of  reaeration  formulae 

Many  authors  use  the  Arrhenius  concept  where 

kt  =  K20eT~™ 

where  0  is  1.024.  Alternatively,  as  in  the  case  of  those  reaeration  formulae  that  use  the  molecular 
diffusivity  of  oxygen,  a  relationship  such  as: 

Dm  (m2  /  s)  =  4.58 E  -117  + 1.2 E  -  9 


where  T  is  the  temperature  in  C,  has  been  used.  These  two  relationships  yield  almost  identical  results 
as  shown  in  Figure  25. 

Upon  further  examination  in  the  current  code  of  W2,  no  temperature  correction  was  made  to  the 
calculated  value  of  KL.  This  will  be  changed  in  the  new  code. 
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Variation  of  KL  with  Temperature 


Figure  25.  Variation  of  KLt/Kl2o  as  a  function  of  temperature. 


Linkage  of  Branches  with  Internal  Head  Boundary 
Conditions 


Linkage  of  mainstem  branches 

One  issue  in  the  development  of  the  river  basin  model  is  the  linkage  of  branches  of  different 
channel  slope  orientation.  This  was  shown  earlier  in  Figure  4.  Figure  26  shows  in  detail  some  of 
the  variable  definitions  with  the  current  sloped  channel  scheme. 
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CE-QUAL-W2  coordinate  system 


p,<E>,P,B 


Figure  26.  Variable  definitions  for  W2  model  with  arbitrary  channel  slope. 


At  an  internal  junction  between  two  branches  where  the  slope  changes,  how  should  the  velocity 
field  and  concentration  field  be  mapped  from  the  upstream  grid  to  the  downstream  grid?  A  first 
order  approach  for  the  velocity  vectors  could  be  to  assume  that  the  angle  orientation  of  the  branch 
does  not  affect  the  stream  lines,  since  they  will  bend  to  the  shape  of  the  channel  bottom  (Figure  27). 
This  is  no  different  than  how  W2  handles  the  transition  between  segments  now  (Figure  28). 


Branch  transition 


Figure  27.  Branch  transition  vertically. 
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longitudinal  velocities  between 
segments  of  varying  orientation 

Figure  28.  Variation  of  longitudinal  velocity  between  segments. 

In  contrast  to  variation  of  the  flow  horizontally,  the  vertical  and  longitudinal  components  of  the 
branch  may  be  important  for  plunging  flows.  This  will  be  deduced  from  numerical  evaluation.  In 
the  case  of  preserving  the  x  and  z  components  of  the  velocity,  the  transition  could  be  illustrated  as 
in  Figure  29. 


Uid-iu=cosaUid 


Wid=sinaUid 


Figure  29.  Illustration  of  branch  transition. 


However,  the  vertical  velocity  of  a  cell  is  not  determined  at  the  side  edge  of  a  segment,  but  at  the 
bottom  of  the  segment.  In  order  for  all  the  volume  to  be  passed  from  one  cell  to  another,  all  the 
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flow  from  the  downstream  segment  (ID)  should  be  transferred  to  the  upstream  segment  (IU).  Since 
the  model  does  not  assume  strong  vertical  accelerations,  we  may  be  forced  to  neglect  the  vertical 
component  of  velocity  at  this  transition  and  assume  that  the  longitudinal  velocity  entering  segment 
IU  is  Uid-  Calculations  and  numerical  testing  will  be  performed  to  verify  this  assumption. 

Another  issue  is  the  linkage  between  branches  when  the  grid  sizes  are  different  between  the 
upstream  grid  and  the  downstream  grid.  Flow  and  mass  will  be  conserved  at  the  linkage  and  will  be 
computed  internally.  This  spatial  averaging  of  the  flow  (and  velocity),  heat  and  mass  to  preserve 
flow  and  constituent  mass  between  branches  is  illustrated  conceptually  in  Figure  30. 


Figure  30.  Transfer  of  mass  and  momentum  between  main  stem  branches  with  unequal  grid 
spacing. 

Linkage  of  tributary  branches 

The  existing  W2  model  assumes  all  tributary  branches  come  in  at  right  angles  to  the  main  channel. 
In  many  cases,  this  is  appropriate.  This  orientation  (Figure  31)  allows  volume  exchange,  but  no 
momentum  exchange  between  branches.  The  CE-QUAL-RTVl  (Environmental  Laboratory,  1995) 
and  the  EPA  DYNHYD  models  (Ambrose,  et  ah,  1988)  also  neglect  momentum  effects  of  lateral 
tributary  inflows.  For  branches  with  arbitrary  channel  orientation  (as  in  Figure  32),  code  changes 
will  be  made  to  allow  momentum,  in  addition  to  volume  (this  is  accounted  for  in  the  free  surface 
equation  as  q),  to  be  exchanged  between  branches. 

In  this  section,  the  linking  of  these  tributary  branches  with  the  main  stem  and  preserving  momentum 
between  them  will  be  discussed. 
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The  tributary  inflow  can  create  shear  stress  along  both  the  longitudinal  axis  of  the  main  stem  branch 
and  along  the  y-axis  of  the  segment.  In  the  current  model,  this  cross-shear  term  is  neglected  and 
does  not  impact  vertical  mixing.  The  only  vertical  mixing  as  a  result  of  cross-shear  is  from  the 
wind  component  in  the  lateral  direction.  For  this  new  formulation,  the  cross-shear  mixing  will  be 
added  to  the  cross-shear  wind  stress  for  the  computation  involving  the  vertical  eddy  viscosity  and 
vertical  diffusivity.  This  involves  determining  the  y  and  x  velocity  components  of  the  entering 
branch  as  shown  in  Figure  33. 


Figure  33.  Schematic  of  branch  connection. 


Longitudinal  momentum 

The  vector  component  of  velocity  in  the  x-direction  of  the  main  channel,  Ux,  can  be  computed  by 
analysis  of  the  channel  orientations.  This  component  in  the  x-direction  would  be:  Ux=Ucosj3  where 
U  is  the  longitudinal  velocity  of  the  tributary  at  segment  ID  for  the  tributary  branch  and  (3  is  the 
difference  in  the  angle  between  the  main  stem  and  tributary  segments  (see  Figure  34). 
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Figure  34.  Schematic  of  computation  of  x  and  y  velocity  components. 


The  conservation  of  momentum  about  a  control  volume,  the  main  stem  segment,  would  result  in  an 
additional  source  of  momentum.  Lai  (1986)  shows  that  the  correction  to  the  x-momentum  equation 
would  be: 


qBux 

where  q  is  the  lateral  inflow  per  unit  length. 

This  arises  from  re-deriving  the  momentum  equations  and  assuming  that  all  the  fluid  (q)  entering 
the  segment  is  moving  at  the  velocity  Ux.  This  correction  to  the  x-momentum  equation  would  be: 


d  UB  dUUB  dWUB 

dt  dx  dz 


gB  sin  a  +  g  cosaB— — 
ox 


dr]  gcosaB  r  dp 


1  gg£xr 

p  dx 


+ 


1  9Btu 
p  dz 


momentum  from  side  tributaries 


Cross-shear  of  tributary  inflow 

The  y-velocity  coming  into  a  reservoir  also  may  contribute  significantly  to  vertical  mixing.  The  y 
component  of  a  tributary  inflow  is  (see  Figure  33):  Uy=Usinp.  Since  there  is  no  y-momentum 
equation,  the  only  mechanism  for  mixing  energy  with  the  present  formulation  of  the  vertical  shear 

stress  is  the  cross-shear  stress  from  the  wind  given  earlier  as  Tlty  =  CD  pa  Wh2  sin(0,  -02) . 
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This  cross-shear  stress  accounts  for  the  shear  stress  and  mixing  that  results  from  wind  blowing 
across  the  y-axis  of  the  segment.  The  lateral  branch  inflow  at  a  velocity,  Uy,  could  be  thought  of  as 
an  additional  component  of  that  stress  under  the  current  context  of  the  turbulence  closure 
approximations. 

Assuming  that  the  water  in  the  y-direction  has  zero  velocity,  the  additional  shear  stress  could  be 
parameterized  as  an  interfacial  shear: 

T  ~nLT12 

^  yt  rib  Mary  —  P  g  ^  y 

where  f  is  an  interfacial  friction  factor.  For  two-layer  flow  systems,  f  has  been  found  to  be  of  order 
0.01.  The  value  of  f  for  this  non-ideal  approach  could  be  determined  by  numerical  computation. 
Hence,  the  value  of  the  cross-shear  term  would  be  increased  by  a  lateral  tributary  inflow.  This  will 
be  evaluated  by  numerical  experiments  computing  the  magnitude  of  the  cross-shear  term  from  wind 
and  from  lateral  inflow.  A  more  robust  theoretical  approach  may  be  needed  to  account  for  this 
increase  in  lateral  shear,  but  that  may  be  necessary  only  if  the  model  includes  the  y-momentum 
equation. 
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Additional  River  Basin  Model  Issues 


Since  W2  will  be  applied  to  river  basins,  additional  factors  may  need  to  be  addressed  in  the  model. 

Some  of  these  potential  factors  may  be: 

•  impact  of  river  meanders  on  transport  (see  Holley  and  Jirka,  1 986) 

•  impact  of  shallow  embayments  -  shallow  embayments  may  be  necessary  for  the 
continuity  balance  but  may  not  really  fully  participate  in  momentum  transport.  This  issue 
is  also  necessary  to  be  addressed  in  reservoir  systems.  In  W2  currently,  one  can  specify  a 
lateral  branch  with  an  angle  90°  to  the  main  stem  to  account  for  this.  Alternatively,  one 
could  follow  the  approach  of  CE-QUAL-RIV 1 ,  where  storage  areas  are  added  to  a 
segment  that  participate  in  water  quality  and  volume  balance,  but  not  in  the  momentum 
balance. 

•  evaluation  of  revising  the  W2  bathymetry  layout  for  river  channels  -  this  would  include 
allowing  the  user  to  define  a  one-layer  cell  that  has  a  changing  cross-sectional  area  with 
depth.  For  example  in  CE-QUAL-RIV  1,  the  user  can  input  a  rectangular,  triangular, 
trapezoidal,  parabolic,  or  ellipsoidal  cross-section.  In  W2,  the  user  would  have  to  input 
several  vertical  layers  to  have  the  same  shape  as  the  CE-QUAL-RIV  1  one-layer  model. 

In  order  to  improve  the  efficiency  of  the  code,  efforts  will  be  made  to  allow  W2  to  also 
have  one-layer  with  an  arbitrary  channel  shape. 

•  impact  of  constrictions  in  the  river  channel  -  in  river  systems,  the  impact  of  bridge  piers 
or  other  obstacles  in  the  fluid  may  be  significant.  The  x-momentum  equation  may  need  to 
reflect  the  additional  frictional  resistance. 
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Summary 


The  following  tasks  will  be  performed  to  implement  the  above  points: 

•  alter  the  x-momentum  equation  for  a  general  channel  slope  within  the  Version  3  code 

•  alter  the  input  of  data  for  a  sloping  channel  by  adding  the  channel  slope  and  explaining 
how  the  channel  geometry  will  be  entered 

•  allow  for  a  branch  to  be  associated  with  a  given  meteorological  file  or  allow  the  option 
for  altitude  dependent  meteorological  terms  to  be  adjusted 

•  allow  for  cell  variable  or  branch  variable  friction  and  wind  sheltering  coefficient 

•  alter  the  code  to  allow  for  internal  head  boundary  conditions  between  a  main  stem  and  a 
tributary  branch  to  transfer  momentum  properly  (both  longitudinal  momentum  and  the 
additional  cross-shear  terms) 

•  alter  the  code  to  allow  a  transition  between  main  stem  branches  when  the  channel  slope 
and  the  cell  layers  are  different 

Once  these  tasks  are  complete,  the  code  will  be  tested  for  test  cases  in  order  to  evaluate  the  impact 

of  the  various  model  changes.  Numerous  vertical  eddy  viscosity  formulations  will  be  used  to 

compare  model  predictions  to  theory. 
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